LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ssterf.f
Go to the documentation of this file.
1*> \brief \b SSTERF
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SSTERF + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssterf.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssterf.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssterf.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SSTERF( N, D, E, INFO )
22*
23* .. Scalar Arguments ..
24* INTEGER INFO, N
25* ..
26* .. Array Arguments ..
27* REAL D( * ), E( * )
28* ..
29*
30*
31*> \par Purpose:
32* =============
33*>
34*> \verbatim
35*>
36*> SSTERF computes all eigenvalues of a symmetric tridiagonal matrix
37*> using the Pal-Walker-Kahan variant of the QL or QR algorithm.
38*> \endverbatim
39*
40* Arguments:
41* ==========
42*
43*> \param[in] N
44*> \verbatim
45*> N is INTEGER
46*> The order of the matrix. N >= 0.
47*> \endverbatim
48*>
49*> \param[in,out] D
50*> \verbatim
51*> D is REAL array, dimension (N)
52*> On entry, the n diagonal elements of the tridiagonal matrix.
53*> On exit, if INFO = 0, the eigenvalues in ascending order.
54*> \endverbatim
55*>
56*> \param[in,out] E
57*> \verbatim
58*> E is REAL array, dimension (N-1)
59*> On entry, the (n-1) subdiagonal elements of the tridiagonal
60*> matrix.
61*> On exit, E has been destroyed.
62*> \endverbatim
63*>
64*> \param[out] INFO
65*> \verbatim
66*> INFO is INTEGER
67*> = 0: successful exit
68*> < 0: if INFO = -i, the i-th argument had an illegal value
69*> > 0: the algorithm failed to find all of the eigenvalues in
70*> a total of 30*N iterations; if INFO = i, then i
71*> elements of E have not converged to zero.
72*> \endverbatim
73*
74* Authors:
75* ========
76*
77*> \author Univ. of Tennessee
78*> \author Univ. of California Berkeley
79*> \author Univ. of Colorado Denver
80*> \author NAG Ltd.
81*
82*> \ingroup sterf
83*
84* =====================================================================
85 SUBROUTINE ssterf( N, D, E, INFO )
86*
87* -- LAPACK computational routine --
88* -- LAPACK is a software package provided by Univ. of Tennessee, --
89* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
90*
91* .. Scalar Arguments ..
92 INTEGER INFO, N
93* ..
94* .. Array Arguments ..
95 REAL D( * ), E( * )
96* ..
97*
98* =====================================================================
99*
100* .. Parameters ..
101 REAL ZERO, ONE, TWO, THREE
102 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
103 $ three = 3.0e0 )
104 INTEGER MAXIT
105 parameter( maxit = 30 )
106* ..
107* .. Local Scalars ..
108 INTEGER I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M,
109 $ NMAXIT
110 REAL ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
111 $ OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
112 $ SIGMA, SSFMAX, SSFMIN
113* ..
114* .. External Functions ..
115 REAL SLAMCH, SLANST, SLAPY2
116 EXTERNAL slamch, slanst, slapy2
117* ..
118* .. External Subroutines ..
119 EXTERNAL slae2, slascl, slasrt, xerbla
120* ..
121* .. Intrinsic Functions ..
122 INTRINSIC abs, sign, sqrt
123* ..
124* .. Executable Statements ..
125*
126* Test the input parameters.
127*
128 info = 0
129*
130* Quick return if possible
131*
132 IF( n.LT.0 ) THEN
133 info = -1
134 CALL xerbla( 'SSTERF', -info )
135 RETURN
136 END IF
137 IF( n.LE.1 )
138 $ RETURN
139*
140* Determine the unit roundoff for this environment.
141*
142 eps = slamch( 'E' )
143 eps2 = eps**2
144 safmin = slamch( 'S' )
145 safmax = one / safmin
146 ssfmax = sqrt( safmax ) / three
147 ssfmin = sqrt( safmin ) / eps2
148*
149* Compute the eigenvalues of the tridiagonal matrix.
150*
151 nmaxit = n*maxit
152 sigma = zero
153 jtot = 0
154*
155* Determine where the matrix splits and choose QL or QR iteration
156* for each block, according to whether top or bottom diagonal
157* element is smaller.
158*
159 l1 = 1
160*
161 10 CONTINUE
162 IF( l1.GT.n )
163 $ GO TO 170
164 IF( l1.GT.1 )
165 $ e( l1-1 ) = zero
166 DO 20 m = l1, n - 1
167 IF( abs( e( m ) ).LE.( sqrt( abs( d( m ) ) )*
168 $ sqrt( abs( d( m+1 ) ) ) )*eps ) THEN
169 e( m ) = zero
170 GO TO 30
171 END IF
172 20 CONTINUE
173 m = n
174*
175 30 CONTINUE
176 l = l1
177 lsv = l
178 lend = m
179 lendsv = lend
180 l1 = m + 1
181 IF( lend.EQ.l )
182 $ GO TO 10
183*
184* Scale submatrix in rows and columns L to LEND
185*
186 anorm = slanst( 'M', lend-l+1, d( l ), e( l ) )
187 iscale = 0
188 IF( anorm.EQ.zero )
189 $ GO TO 10
190 IF( anorm.GT.ssfmax ) THEN
191 iscale = 1
192 CALL slascl( 'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
193 $ info )
194 CALL slascl( 'G', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
195 $ info )
196 ELSE IF( anorm.LT.ssfmin ) THEN
197 iscale = 2
198 CALL slascl( 'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
199 $ info )
200 CALL slascl( 'G', 0, 0, anorm, ssfmin, lend-l, 1, e( l ), n,
201 $ info )
202 END IF
203*
204 DO 40 i = l, lend - 1
205 e( i ) = e( i )**2
206 40 CONTINUE
207*
208* Choose between QL and QR iteration
209*
210 IF( abs( d( lend ) ).LT.abs( d( l ) ) ) THEN
211 lend = lsv
212 l = lendsv
213 END IF
214*
215 IF( lend.GE.l ) THEN
216*
217* QL Iteration
218*
219* Look for small subdiagonal element.
220*
221 50 CONTINUE
222 IF( l.NE.lend ) THEN
223 DO 60 m = l, lend - 1
224 IF( abs( e( m ) ).LE.eps2*abs( d( m )*d( m+1 ) ) )
225 $ GO TO 70
226 60 CONTINUE
227 END IF
228 m = lend
229*
230 70 CONTINUE
231 IF( m.LT.lend )
232 $ e( m ) = zero
233 p = d( l )
234 IF( m.EQ.l )
235 $ GO TO 90
236*
237* If remaining matrix is 2 by 2, use SLAE2 to compute its
238* eigenvalues.
239*
240 IF( m.EQ.l+1 ) THEN
241 rte = sqrt( e( l ) )
242 CALL slae2( d( l ), rte, d( l+1 ), rt1, rt2 )
243 d( l ) = rt1
244 d( l+1 ) = rt2
245 e( l ) = zero
246 l = l + 2
247 IF( l.LE.lend )
248 $ GO TO 50
249 GO TO 150
250 END IF
251*
252 IF( jtot.EQ.nmaxit )
253 $ GO TO 150
254 jtot = jtot + 1
255*
256* Form shift.
257*
258 rte = sqrt( e( l ) )
259 sigma = ( d( l+1 )-p ) / ( two*rte )
260 r = slapy2( sigma, one )
261 sigma = p - ( rte / ( sigma+sign( r, sigma ) ) )
262*
263 c = one
264 s = zero
265 gamma = d( m ) - sigma
266 p = gamma*gamma
267*
268* Inner loop
269*
270 DO 80 i = m - 1, l, -1
271 bb = e( i )
272 r = p + bb
273 IF( i.NE.m-1 )
274 $ e( i+1 ) = s*r
275 oldc = c
276 c = p / r
277 s = bb / r
278 oldgam = gamma
279 alpha = d( i )
280 gamma = c*( alpha-sigma ) - s*oldgam
281 d( i+1 ) = oldgam + ( alpha-gamma )
282 IF( c.NE.zero ) THEN
283 p = ( gamma*gamma ) / c
284 ELSE
285 p = oldc*bb
286 END IF
287 80 CONTINUE
288*
289 e( l ) = s*p
290 d( l ) = sigma + gamma
291 GO TO 50
292*
293* Eigenvalue found.
294*
295 90 CONTINUE
296 d( l ) = p
297*
298 l = l + 1
299 IF( l.LE.lend )
300 $ GO TO 50
301 GO TO 150
302*
303 ELSE
304*
305* QR Iteration
306*
307* Look for small superdiagonal element.
308*
309 100 CONTINUE
310 DO 110 m = l, lend + 1, -1
311 IF( abs( e( m-1 ) ).LE.eps2*abs( d( m )*d( m-1 ) ) )
312 $ GO TO 120
313 110 CONTINUE
314 m = lend
315*
316 120 CONTINUE
317 IF( m.GT.lend )
318 $ e( m-1 ) = zero
319 p = d( l )
320 IF( m.EQ.l )
321 $ GO TO 140
322*
323* If remaining matrix is 2 by 2, use SLAE2 to compute its
324* eigenvalues.
325*
326 IF( m.EQ.l-1 ) THEN
327 rte = sqrt( e( l-1 ) )
328 CALL slae2( d( l ), rte, d( l-1 ), rt1, rt2 )
329 d( l ) = rt1
330 d( l-1 ) = rt2
331 e( l-1 ) = zero
332 l = l - 2
333 IF( l.GE.lend )
334 $ GO TO 100
335 GO TO 150
336 END IF
337*
338 IF( jtot.EQ.nmaxit )
339 $ GO TO 150
340 jtot = jtot + 1
341*
342* Form shift.
343*
344 rte = sqrt( e( l-1 ) )
345 sigma = ( d( l-1 )-p ) / ( two*rte )
346 r = slapy2( sigma, one )
347 sigma = p - ( rte / ( sigma+sign( r, sigma ) ) )
348*
349 c = one
350 s = zero
351 gamma = d( m ) - sigma
352 p = gamma*gamma
353*
354* Inner loop
355*
356 DO 130 i = m, l - 1
357 bb = e( i )
358 r = p + bb
359 IF( i.NE.m )
360 $ e( i-1 ) = s*r
361 oldc = c
362 c = p / r
363 s = bb / r
364 oldgam = gamma
365 alpha = d( i+1 )
366 gamma = c*( alpha-sigma ) - s*oldgam
367 d( i ) = oldgam + ( alpha-gamma )
368 IF( c.NE.zero ) THEN
369 p = ( gamma*gamma ) / c
370 ELSE
371 p = oldc*bb
372 END IF
373 130 CONTINUE
374*
375 e( l-1 ) = s*p
376 d( l ) = sigma + gamma
377 GO TO 100
378*
379* Eigenvalue found.
380*
381 140 CONTINUE
382 d( l ) = p
383*
384 l = l - 1
385 IF( l.GE.lend )
386 $ GO TO 100
387 GO TO 150
388*
389 END IF
390*
391* Undo scaling if necessary
392*
393 150 CONTINUE
394 IF( iscale.EQ.1 )
395 $ CALL slascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
396 $ d( lsv ), n, info )
397 IF( iscale.EQ.2 )
398 $ CALL slascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
399 $ d( lsv ), n, info )
400*
401* Check for no convergence to an eigenvalue after a total
402* of N*MAXIT iterations.
403*
404 IF( jtot.LT.nmaxit )
405 $ GO TO 10
406 DO 160 i = 1, n - 1
407 IF( e( i ).NE.zero )
408 $ info = info + 1
409 160 CONTINUE
410 GO TO 180
411*
412* Sort eigenvalues in increasing order.
413*
414 170 CONTINUE
415 CALL slasrt( 'I', n, d, info )
416*
417 180 CONTINUE
418 RETURN
419*
420* End of SSTERF
421*
422 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
slae2
subroutine slae2(a, b, c, rt1, rt2)
SLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition slae2.f:102
slascl
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition slascl.f:143
slasrt
subroutine slasrt(id, n, d, info)
SLASRT sorts numbers in increasing or decreasing order.
Definition slasrt.f:88
ssterf
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:86