LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 *> \ingroup auxOTHERcomputational
83 *
84 * =====================================================================
85  SUBROUTINE dsterf( 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  DOUBLE PRECISION D( * ), E( * )
96 * ..
97 *
98 * =====================================================================
99 *
100 * .. Parameters ..
101  DOUBLE PRECISION ZERO, ONE, TWO, THREE
102  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
103  $ three = 3.0d0 )
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  DOUBLE PRECISION ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
111  $ OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
112  $ SIGMA, SSFMAX, SSFMIN, RMAX
113 * ..
114 * .. External Functions ..
115  DOUBLE PRECISION DLAMCH, DLANST, DLAPY2
116  EXTERNAL dlamch, dlanst, dlapy2
117 * ..
118 * .. External Subroutines ..
119  EXTERNAL dlae2, dlascl, dlasrt, 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( 'DSTERF', -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 = dlamch( 'E' )
143  eps2 = eps**2
144  safmin = dlamch( 'S' )
145  safmax = one / safmin
146  ssfmax = sqrt( safmax ) / three
147  ssfmin = sqrt( safmin ) / eps2
148  rmax = dlamch( 'O' )
149 *
150 * Compute the eigenvalues of the tridiagonal matrix.
151 *
152  nmaxit = n*maxit
153  sigma = zero
154  jtot = 0
155 *
156 * Determine where the matrix splits and choose QL or QR iteration
157 * for each block, according to whether top or bottom diagonal
158 * element is smaller.
159 *
160  l1 = 1
161 *
162  10 CONTINUE
163  IF( l1.GT.n )
164  $ GO TO 170
165  IF( l1.GT.1 )
166  $ e( l1-1 ) = zero
167  DO 20 m = l1, n - 1
168  IF( abs( e( m ) ).LE.( sqrt( abs( d( m ) ) )*sqrt( abs( d( m+
169  $ 1 ) ) ) )*eps ) THEN
170  e( m ) = zero
171  GO TO 30
172  END IF
173  20 CONTINUE
174  m = n
175 *
176  30 CONTINUE
177  l = l1
178  lsv = l
179  lend = m
180  lendsv = lend
181  l1 = m + 1
182  IF( lend.EQ.l )
183  $ GO TO 10
184 *
185 * Scale submatrix in rows and columns L to LEND
186 *
187  anorm = dlanst( 'M', lend-l+1, d( l ), e( l ) )
188  iscale = 0
189  IF( anorm.EQ.zero )
190  $ GO TO 10
191  IF( (anorm.GT.ssfmax) ) THEN
192  iscale = 1
193  CALL dlascl( 'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
194  $ info )
195  CALL dlascl( 'G', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
196  $ info )
197  ELSE IF( anorm.LT.ssfmin ) THEN
198  iscale = 2
199  CALL dlascl( 'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
200  $ info )
201  CALL dlascl( 'G', 0, 0, anorm, ssfmin, lend-l, 1, e( l ), n,
202  $ info )
203  END IF
204 *
205  DO 40 i = l, lend - 1
206  e( i ) = e( i )**2
207  40 CONTINUE
208 *
209 * Choose between QL and QR iteration
210 *
211  IF( abs( d( lend ) ).LT.abs( d( l ) ) ) THEN
212  lend = lsv
213  l = lendsv
214  END IF
215 *
216  IF( lend.GE.l ) THEN
217 *
218 * QL Iteration
219 *
220 * Look for small subdiagonal element.
221 *
222  50 CONTINUE
223  IF( l.NE.lend ) THEN
224  DO 60 m = l, lend - 1
225  IF( abs( e( m ) ).LE.eps2*abs( d( m )*d( m+1 ) ) )
226  $ GO TO 70
227  60 CONTINUE
228  END IF
229  m = lend
230 *
231  70 CONTINUE
232  IF( m.LT.lend )
233  $ e( m ) = zero
234  p = d( l )
235  IF( m.EQ.l )
236  $ GO TO 90
237 *
238 * If remaining matrix is 2 by 2, use DLAE2 to compute its
239 * eigenvalues.
240 *
241  IF( m.EQ.l+1 ) THEN
242  rte = sqrt( e( l ) )
243  CALL dlae2( d( l ), rte, d( l+1 ), rt1, rt2 )
244  d( l ) = rt1
245  d( l+1 ) = rt2
246  e( l ) = zero
247  l = l + 2
248  IF( l.LE.lend )
249  $ GO TO 50
250  GO TO 150
251  END IF
252 *
253  IF( jtot.EQ.nmaxit )
254  $ GO TO 150
255  jtot = jtot + 1
256 *
257 * Form shift.
258 *
259  rte = sqrt( e( l ) )
260  sigma = ( d( l+1 )-p ) / ( two*rte )
261  r = dlapy2( sigma, one )
262  sigma = p - ( rte / ( sigma+sign( r, sigma ) ) )
263 *
264  c = one
265  s = zero
266  gamma = d( m ) - sigma
267  p = gamma*gamma
268 *
269 * Inner loop
270 *
271  DO 80 i = m - 1, l, -1
272  bb = e( i )
273  r = p + bb
274  IF( i.NE.m-1 )
275  $ e( i+1 ) = s*r
276  oldc = c
277  c = p / r
278  s = bb / r
279  oldgam = gamma
280  alpha = d( i )
281  gamma = c*( alpha-sigma ) - s*oldgam
282  d( i+1 ) = oldgam + ( alpha-gamma )
283  IF( c.NE.zero ) THEN
284  p = ( gamma*gamma ) / c
285  ELSE
286  p = oldc*bb
287  END IF
288  80 CONTINUE
289 *
290  e( l ) = s*p
291  d( l ) = sigma + gamma
292  GO TO 50
293 *
294 * Eigenvalue found.
295 *
296  90 CONTINUE
297  d( l ) = p
298 *
299  l = l + 1
300  IF( l.LE.lend )
301  $ GO TO 50
302  GO TO 150
303 *
304  ELSE
305 *
306 * QR Iteration
307 *
308 * Look for small superdiagonal element.
309 *
310  100 CONTINUE
311  DO 110 m = l, lend + 1, -1
312  IF( abs( e( m-1 ) ).LE.eps2*abs( d( m )*d( m-1 ) ) )
313  $ GO TO 120
314  110 CONTINUE
315  m = lend
316 *
317  120 CONTINUE
318  IF( m.GT.lend )
319  $ e( m-1 ) = zero
320  p = d( l )
321  IF( m.EQ.l )
322  $ GO TO 140
323 *
324 * If remaining matrix is 2 by 2, use DLAE2 to compute its
325 * eigenvalues.
326 *
327  IF( m.EQ.l-1 ) THEN
328  rte = sqrt( e( l-1 ) )
329  CALL dlae2( d( l ), rte, d( l-1 ), rt1, rt2 )
330  d( l ) = rt1
331  d( l-1 ) = rt2
332  e( l-1 ) = zero
333  l = l - 2
334  IF( l.GE.lend )
335  $ GO TO 100
336  GO TO 150
337  END IF
338 *
339  IF( jtot.EQ.nmaxit )
340  $ GO TO 150
341  jtot = jtot + 1
342 *
343 * Form shift.
344 *
345  rte = sqrt( e( l-1 ) )
346  sigma = ( d( l-1 )-p ) / ( two*rte )
347  r = dlapy2( sigma, one )
348  sigma = p - ( rte / ( sigma+sign( r, sigma ) ) )
349 *
350  c = one
351  s = zero
352  gamma = d( m ) - sigma
353  p = gamma*gamma
354 *
355 * Inner loop
356 *
357  DO 130 i = m, l - 1
358  bb = e( i )
359  r = p + bb
360  IF( i.NE.m )
361  $ e( i-1 ) = s*r
362  oldc = c
363  c = p / r
364  s = bb / r
365  oldgam = gamma
366  alpha = d( i+1 )
367  gamma = c*( alpha-sigma ) - s*oldgam
368  d( i ) = oldgam + ( alpha-gamma )
369  IF( c.NE.zero ) THEN
370  p = ( gamma*gamma ) / c
371  ELSE
372  p = oldc*bb
373  END IF
374  130 CONTINUE
375 *
376  e( l-1 ) = s*p
377  d( l ) = sigma + gamma
378  GO TO 100
379 *
380 * Eigenvalue found.
381 *
382  140 CONTINUE
383  d( l ) = p
384 *
385  l = l - 1
386  IF( l.GE.lend )
387  $ GO TO 100
388  GO TO 150
389 *
390  END IF
391 *
392 * Undo scaling if necessary
393 *
394  150 CONTINUE
395  IF( iscale.EQ.1 )
396  $ CALL dlascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
397  $ d( lsv ), n, info )
398  IF( iscale.EQ.2 )
399  $ CALL dlascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
400  $ d( lsv ), n, info )
401 *
402 * Check for no convergence to an eigenvalue after a total
403 * of N*MAXIT iterations.
404 *
405  IF( jtot.LT.nmaxit )
406  $ GO TO 10
407  DO 160 i = 1, n - 1
408  IF( e( i ).NE.zero )
409  $ info = info + 1
410  160 CONTINUE
411  GO TO 180
412 *
413 * Sort eigenvalues in increasing order.
414 *
415  170 CONTINUE
416  CALL dlasrt( 'I', n, d, info )
417 *
418  180 CONTINUE
419  RETURN
420 *
421 * End of DSTERF
422 *
423  END
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:143
subroutine dlae2(A, B, C, RT1, RT2)
DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition: dlae2.f:102
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:88
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86