LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dsterf.f
Go to the documentation of this file.
1 *> \brief \b DSTERF
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DSTERF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsterf.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsterf.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsterf.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DSTERF( N, D, E, INFO )
22 *
23 * .. Scalar Arguments ..
24 * INTEGER INFO, N
25 * ..
26 * .. Array Arguments ..
27 * DOUBLE PRECISION D( * ), E( * )
28 * ..
29 *
30 *
31 *> \par Purpose:
32 * =============
33 *>
34 *> \verbatim
35 *>
36 *> DSTERF 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 dsterf( 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  DOUBLE PRECISION d( * ), e( * )
99 * ..
100 *
101 * =====================================================================
102 *
103 * .. Parameters ..
104  DOUBLE PRECISION zero, one, two, three
105  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
106  $ three = 3.0d0 )
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  DOUBLE PRECISION alpha, anorm, bb, c, eps, eps2, gamma, oldc,
114  $ oldgam, p, r, rt1, rt2, rte, s, safmax, safmin,
115  $ sigma, ssfmax, ssfmin, rmax
116 * ..
117 * .. External Functions ..
118  DOUBLE PRECISION dlamch, dlanst, dlapy2
119  EXTERNAL dlamch, dlanst, dlapy2
120 * ..
121 * .. External Subroutines ..
122  EXTERNAL dlae2, dlascl, dlasrt, 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( 'DSTERF', -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 = dlamch( 'E' )
146  eps2 = eps**2
147  safmin = dlamch( 'S' )
148  safmax = one / safmin
149  ssfmax = sqrt( safmax ) / three
150  ssfmin = sqrt( safmin ) / eps2
151  rmax = dlamch( 'O' )
152 *
153 * Compute the eigenvalues of the tridiagonal matrix.
154 *
155  nmaxit = n*maxit
156  sigma = zero
157  jtot = 0
158 *
159 * Determine where the matrix splits and choose QL or QR iteration
160 * for each block, according to whether top or bottom diagonal
161 * element is smaller.
162 *
163  l1 = 1
164 *
165  10 continue
166  IF( l1.GT.n )
167  $ go to 170
168  IF( l1.GT.1 )
169  $ e( l1-1 ) = zero
170  DO 20 m = l1, n - 1
171  IF( abs( e( m ) ).LE.( sqrt( abs( d( m ) ) )*sqrt( abs( d( m+
172  $ 1 ) ) ) )*eps ) THEN
173  e( m ) = zero
174  go to 30
175  END IF
176  20 continue
177  m = n
178 *
179  30 continue
180  l = l1
181  lsv = l
182  lend = m
183  lendsv = lend
184  l1 = m + 1
185  IF( lend.EQ.l )
186  $ go to 10
187 *
188 * Scale submatrix in rows and columns L to LEND
189 *
190  anorm = dlanst( 'M', lend-l+1, d( l ), e( l ) )
191  iscale = 0
192  IF( anorm.EQ.zero )
193  $ go to 10
194  IF( (anorm.GT.ssfmax) ) THEN
195  iscale = 1
196  CALL dlascl( 'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
197  $ info )
198  CALL dlascl( 'G', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
199  $ info )
200  ELSE IF( anorm.LT.ssfmin ) THEN
201  iscale = 2
202  CALL dlascl( 'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
203  $ info )
204  CALL dlascl( 'G', 0, 0, anorm, ssfmin, lend-l, 1, e( l ), n,
205  $ info )
206  END IF
207 *
208  DO 40 i = l, lend - 1
209  e( i ) = e( i )**2
210  40 continue
211 *
212 * Choose between QL and QR iteration
213 *
214  IF( abs( d( lend ) ).LT.abs( d( l ) ) ) THEN
215  lend = lsv
216  l = lendsv
217  END IF
218 *
219  IF( lend.GE.l ) THEN
220 *
221 * QL Iteration
222 *
223 * Look for small subdiagonal element.
224 *
225  50 continue
226  IF( l.NE.lend ) THEN
227  DO 60 m = l, lend - 1
228  IF( abs( e( m ) ).LE.eps2*abs( d( m )*d( m+1 ) ) )
229  $ go to 70
230  60 continue
231  END IF
232  m = lend
233 *
234  70 continue
235  IF( m.LT.lend )
236  $ e( m ) = zero
237  p = d( l )
238  IF( m.EQ.l )
239  $ go to 90
240 *
241 * If remaining matrix is 2 by 2, use DLAE2 to compute its
242 * eigenvalues.
243 *
244  IF( m.EQ.l+1 ) THEN
245  rte = sqrt( e( l ) )
246  CALL dlae2( d( l ), rte, d( l+1 ), rt1, rt2 )
247  d( l ) = rt1
248  d( l+1 ) = rt2
249  e( l ) = zero
250  l = l + 2
251  IF( l.LE.lend )
252  $ go to 50
253  go to 150
254  END IF
255 *
256  IF( jtot.EQ.nmaxit )
257  $ go to 150
258  jtot = jtot + 1
259 *
260 * Form shift.
261 *
262  rte = sqrt( e( l ) )
263  sigma = ( d( l+1 )-p ) / ( two*rte )
264  r = dlapy2( sigma, one )
265  sigma = p - ( rte / ( sigma+sign( r, sigma ) ) )
266 *
267  c = one
268  s = zero
269  gamma = d( m ) - sigma
270  p = gamma*gamma
271 *
272 * Inner loop
273 *
274  DO 80 i = m - 1, l, -1
275  bb = e( i )
276  r = p + bb
277  IF( i.NE.m-1 )
278  $ e( i+1 ) = s*r
279  oldc = c
280  c = p / r
281  s = bb / r
282  oldgam = gamma
283  alpha = d( i )
284  gamma = c*( alpha-sigma ) - s*oldgam
285  d( i+1 ) = oldgam + ( alpha-gamma )
286  IF( c.NE.zero ) THEN
287  p = ( gamma*gamma ) / c
288  ELSE
289  p = oldc*bb
290  END IF
291  80 continue
292 *
293  e( l ) = s*p
294  d( l ) = sigma + gamma
295  go to 50
296 *
297 * Eigenvalue found.
298 *
299  90 continue
300  d( l ) = p
301 *
302  l = l + 1
303  IF( l.LE.lend )
304  $ go to 50
305  go to 150
306 *
307  ELSE
308 *
309 * QR Iteration
310 *
311 * Look for small superdiagonal element.
312 *
313  100 continue
314  DO 110 m = l, lend + 1, -1
315  IF( abs( e( m-1 ) ).LE.eps2*abs( d( m )*d( m-1 ) ) )
316  $ go to 120
317  110 continue
318  m = lend
319 *
320  120 continue
321  IF( m.GT.lend )
322  $ e( m-1 ) = zero
323  p = d( l )
324  IF( m.EQ.l )
325  $ go to 140
326 *
327 * If remaining matrix is 2 by 2, use DLAE2 to compute its
328 * eigenvalues.
329 *
330  IF( m.EQ.l-1 ) THEN
331  rte = sqrt( e( l-1 ) )
332  CALL dlae2( d( l ), rte, d( l-1 ), rt1, rt2 )
333  d( l ) = rt1
334  d( l-1 ) = rt2
335  e( l-1 ) = zero
336  l = l - 2
337  IF( l.GE.lend )
338  $ go to 100
339  go to 150
340  END IF
341 *
342  IF( jtot.EQ.nmaxit )
343  $ go to 150
344  jtot = jtot + 1
345 *
346 * Form shift.
347 *
348  rte = sqrt( e( l-1 ) )
349  sigma = ( d( l-1 )-p ) / ( two*rte )
350  r = dlapy2( sigma, one )
351  sigma = p - ( rte / ( sigma+sign( r, sigma ) ) )
352 *
353  c = one
354  s = zero
355  gamma = d( m ) - sigma
356  p = gamma*gamma
357 *
358 * Inner loop
359 *
360  DO 130 i = m, l - 1
361  bb = e( i )
362  r = p + bb
363  IF( i.NE.m )
364  $ e( i-1 ) = s*r
365  oldc = c
366  c = p / r
367  s = bb / r
368  oldgam = gamma
369  alpha = d( i+1 )
370  gamma = c*( alpha-sigma ) - s*oldgam
371  d( i ) = oldgam + ( alpha-gamma )
372  IF( c.NE.zero ) THEN
373  p = ( gamma*gamma ) / c
374  ELSE
375  p = oldc*bb
376  END IF
377  130 continue
378 *
379  e( l-1 ) = s*p
380  d( l ) = sigma + gamma
381  go to 100
382 *
383 * Eigenvalue found.
384 *
385  140 continue
386  d( l ) = p
387 *
388  l = l - 1
389  IF( l.GE.lend )
390  $ go to 100
391  go to 150
392 *
393  END IF
394 *
395 * Undo scaling if necessary
396 *
397  150 continue
398  IF( iscale.EQ.1 )
399  $ CALL dlascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
400  $ d( lsv ), n, info )
401  IF( iscale.EQ.2 )
402  $ CALL dlascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
403  $ d( lsv ), n, info )
404 *
405 * Check for no convergence to an eigenvalue after a total
406 * of N*MAXIT iterations.
407 *
408  IF( jtot.LT.nmaxit )
409  $ go to 10
410  DO 160 i = 1, n - 1
411  IF( e( i ).NE.zero )
412  $ info = info + 1
413  160 continue
414  go to 180
415 *
416 * Sort eigenvalues in increasing order.
417 *
418  170 continue
419  CALL dlasrt( 'I', n, d, info )
420 *
421  180 continue
422  return
423 *
424 * End of DSTERF
425 *
426  END