LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date November 2011
83 *
84 *> \ingroup auxOTHERcomputational
85 *
86 * =====================================================================
87  SUBROUTINE ssterf( N, D, E, INFO )
88 *
89 * -- LAPACK computational routine (version 3.4.0) --
90 * -- LAPACK is a software package provided by Univ. of Tennessee, --
91 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
92 * November 2011
93 *
94 * .. Scalar Arguments ..
95  INTEGER info, n
96 * ..
97 * .. Array Arguments ..
98  REAL d( * ), e( * )
99 * ..
100 *
101 * =====================================================================
102 *
103 * .. Parameters ..
104  REAL zero, one, two, three
105  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
106  $ three = 3.0e0 )
107  INTEGER maxit
108  parameter( maxit = 30 )
109 * ..
110 * .. Local Scalars ..
111  INTEGER i, iscale, jtot, l, l1, lend, lendsv, lsv, m,
112  $ nmaxit
113  REAL alpha, anorm, bb, c, eps, eps2, gamma, oldc,
114  $ oldgam, p, r, rt1, rt2, rte, s, safmax, safmin,
115  $ sigma, ssfmax, ssfmin
116 * ..
117 * .. External Functions ..
118  REAL slamch, slanst, slapy2
119  EXTERNAL slamch, slanst, slapy2
120 * ..
121 * .. External Subroutines ..
122  EXTERNAL slae2, slascl, slasrt, xerbla
123 * ..
124 * .. Intrinsic Functions ..
125  INTRINSIC abs, sign, sqrt
126 * ..
127 * .. Executable Statements ..
128 *
129 * Test the input parameters.
130 *
131  info = 0
132 *
133 * Quick return if possible
134 *
135  IF( n.LT.0 ) THEN
136  info = -1
137  CALL xerbla( 'SSTERF', -info )
138  return
139  END IF
140  IF( n.LE.1 )
141  $ return
142 *
143 * Determine the unit roundoff for this environment.
144 *
145  eps = slamch( 'E' )
146  eps2 = eps**2
147  safmin = slamch( 'S' )
148  safmax = one / safmin
149  ssfmax = sqrt( safmax ) / three
150  ssfmin = sqrt( safmin ) / eps2
151 *
152 * Compute the eigenvalues of the tridiagonal matrix.
153 *
154  nmaxit = n*maxit
155  sigma = zero
156  jtot = 0
157 *
158 * Determine where the matrix splits and choose QL or QR iteration
159 * for each block, according to whether top or bottom diagonal
160 * element is smaller.
161 *
162  l1 = 1
163 *
164  10 continue
165  IF( l1.GT.n )
166  $ go to 170
167  IF( l1.GT.1 )
168  $ e( l1-1 ) = zero
169  DO 20 m = l1, n - 1
170  IF( abs( e( m ) ).LE.( sqrt( abs( d( m ) ) )*
171  $ sqrt( abs( d( m+1 ) ) ) )*eps ) THEN
172  e( m ) = zero
173  go to 30
174  END IF
175  20 continue
176  m = n
177 *
178  30 continue
179  l = l1
180  lsv = l
181  lend = m
182  lendsv = lend
183  l1 = m + 1
184  IF( lend.EQ.l )
185  $ go to 10
186 *
187 * Scale submatrix in rows and columns L to LEND
188 *
189  anorm = slanst( 'M', lend-l+1, d( l ), e( l ) )
190  iscale = 0
191  IF( anorm.EQ.zero )
192  $ go to 10
193  IF( anorm.GT.ssfmax ) THEN
194  iscale = 1
195  CALL slascl( 'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
196  $ info )
197  CALL slascl( 'G', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
198  $ info )
199  ELSE IF( anorm.LT.ssfmin ) THEN
200  iscale = 2
201  CALL slascl( 'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
202  $ info )
203  CALL slascl( 'G', 0, 0, anorm, ssfmin, lend-l, 1, e( l ), n,
204  $ info )
205  END IF
206 *
207  DO 40 i = l, lend - 1
208  e( i ) = e( i )**2
209  40 continue
210 *
211 * Choose between QL and QR iteration
212 *
213  IF( abs( d( lend ) ).LT.abs( d( l ) ) ) THEN
214  lend = lsv
215  l = lendsv
216  END IF
217 *
218  IF( lend.GE.l ) THEN
219 *
220 * QL Iteration
221 *
222 * Look for small subdiagonal element.
223 *
224  50 continue
225  IF( l.NE.lend ) THEN
226  DO 60 m = l, lend - 1
227  IF( abs( e( m ) ).LE.eps2*abs( d( m )*d( m+1 ) ) )
228  $ go to 70
229  60 continue
230  END IF
231  m = lend
232 *
233  70 continue
234  IF( m.LT.lend )
235  $ e( m ) = zero
236  p = d( l )
237  IF( m.EQ.l )
238  $ go to 90
239 *
240 * If remaining matrix is 2 by 2, use SLAE2 to compute its
241 * eigenvalues.
242 *
243  IF( m.EQ.l+1 ) THEN
244  rte = sqrt( e( l ) )
245  CALL slae2( d( l ), rte, d( l+1 ), rt1, rt2 )
246  d( l ) = rt1
247  d( l+1 ) = rt2
248  e( l ) = zero
249  l = l + 2
250  IF( l.LE.lend )
251  $ go to 50
252  go to 150
253  END IF
254 *
255  IF( jtot.EQ.nmaxit )
256  $ go to 150
257  jtot = jtot + 1
258 *
259 * Form shift.
260 *
261  rte = sqrt( e( l ) )
262  sigma = ( d( l+1 )-p ) / ( two*rte )
263  r = slapy2( sigma, one )
264  sigma = p - ( rte / ( sigma+sign( r, sigma ) ) )
265 *
266  c = one
267  s = zero
268  gamma = d( m ) - sigma
269  p = gamma*gamma
270 *
271 * Inner loop
272 *
273  DO 80 i = m - 1, l, -1
274  bb = e( i )
275  r = p + bb
276  IF( i.NE.m-1 )
277  $ e( i+1 ) = s*r
278  oldc = c
279  c = p / r
280  s = bb / r
281  oldgam = gamma
282  alpha = d( i )
283  gamma = c*( alpha-sigma ) - s*oldgam
284  d( i+1 ) = oldgam + ( alpha-gamma )
285  IF( c.NE.zero ) THEN
286  p = ( gamma*gamma ) / c
287  ELSE
288  p = oldc*bb
289  END IF
290  80 continue
291 *
292  e( l ) = s*p
293  d( l ) = sigma + gamma
294  go to 50
295 *
296 * Eigenvalue found.
297 *
298  90 continue
299  d( l ) = p
300 *
301  l = l + 1
302  IF( l.LE.lend )
303  $ go to 50
304  go to 150
305 *
306  ELSE
307 *
308 * QR Iteration
309 *
310 * Look for small superdiagonal element.
311 *
312  100 continue
313  DO 110 m = l, lend + 1, -1
314  IF( abs( e( m-1 ) ).LE.eps2*abs( d( m )*d( m-1 ) ) )
315  $ go to 120
316  110 continue
317  m = lend
318 *
319  120 continue
320  IF( m.GT.lend )
321  $ e( m-1 ) = zero
322  p = d( l )
323  IF( m.EQ.l )
324  $ go to 140
325 *
326 * If remaining matrix is 2 by 2, use SLAE2 to compute its
327 * eigenvalues.
328 *
329  IF( m.EQ.l-1 ) THEN
330  rte = sqrt( e( l-1 ) )
331  CALL slae2( d( l ), rte, d( l-1 ), rt1, rt2 )
332  d( l ) = rt1
333  d( l-1 ) = rt2
334  e( l-1 ) = zero
335  l = l - 2
336  IF( l.GE.lend )
337  $ go to 100
338  go to 150
339  END IF
340 *
341  IF( jtot.EQ.nmaxit )
342  $ go to 150
343  jtot = jtot + 1
344 *
345 * Form shift.
346 *
347  rte = sqrt( e( l-1 ) )
348  sigma = ( d( l-1 )-p ) / ( two*rte )
349  r = slapy2( sigma, one )
350  sigma = p - ( rte / ( sigma+sign( r, sigma ) ) )
351 *
352  c = one
353  s = zero
354  gamma = d( m ) - sigma
355  p = gamma*gamma
356 *
357 * Inner loop
358 *
359  DO 130 i = m, l - 1
360  bb = e( i )
361  r = p + bb
362  IF( i.NE.m )
363  $ e( i-1 ) = s*r
364  oldc = c
365  c = p / r
366  s = bb / r
367  oldgam = gamma
368  alpha = d( i+1 )
369  gamma = c*( alpha-sigma ) - s*oldgam
370  d( i ) = oldgam + ( alpha-gamma )
371  IF( c.NE.zero ) THEN
372  p = ( gamma*gamma ) / c
373  ELSE
374  p = oldc*bb
375  END IF
376  130 continue
377 *
378  e( l-1 ) = s*p
379  d( l ) = sigma + gamma
380  go to 100
381 *
382 * Eigenvalue found.
383 *
384  140 continue
385  d( l ) = p
386 *
387  l = l - 1
388  IF( l.GE.lend )
389  $ go to 100
390  go to 150
391 *
392  END IF
393 *
394 * Undo scaling if necessary
395 *
396  150 continue
397  IF( iscale.EQ.1 )
398  $ CALL slascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
399  $ d( lsv ), n, info )
400  IF( iscale.EQ.2 )
401  $ CALL slascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
402  $ d( lsv ), n, info )
403 *
404 * Check for no convergence to an eigenvalue after a total
405 * of N*MAXIT iterations.
406 *
407  IF( jtot.LT.nmaxit )
408  $ go to 10
409  DO 160 i = 1, n - 1
410  IF( e( i ).NE.zero )
411  $ info = info + 1
412  160 continue
413  go to 180
414 *
415 * Sort eigenvalues in increasing order.
416 *
417  170 continue
418  CALL slasrt( 'I', n, d, info )
419 *
420  180 continue
421  return
422 *
423 * End of SSTERF
424 *
425  END