LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dsterf ( integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
integer  INFO 
)

DSTERF

Download DSTERF + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
 DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
 using the Pal-Walker-Kahan variant of the QL or QR algorithm.
Parameters
[in]N
          N is INTEGER
          The order of the matrix.  N >= 0.
[in,out]D
          D is DOUBLE PRECISION array, dimension (N)
          On entry, the n diagonal elements of the tridiagonal matrix.
          On exit, if INFO = 0, the eigenvalues in ascending order.
[in,out]E
          E is DOUBLE PRECISION array, dimension (N-1)
          On entry, the (n-1) subdiagonal elements of the tridiagonal
          matrix.
          On exit, E has been destroyed.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  the algorithm failed to find all of the eigenvalues in
                a total of 30*N iterations; if INFO = i, then i
                elements of E have not converged to zero.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 88 of file dsterf.f.

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 *
dlasrt
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:90
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
dlae2
subroutine dlae2(A, B, C, RT1, RT2)
DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition: dlae2.f:104
dlascl
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:145
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dlapy2
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:65
dlanst
double precision function dlanst(NORM, N, D, E)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix.
Definition: dlanst.f:102

Here is the call graph for this function:

Here is the caller graph for this function: