LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ ssterf()

subroutine ssterf ( integer  N,
real, dimension( * )  D,
real, dimension( * )  E,
integer  INFO 
)

SSTERF

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

Purpose:
 SSTERF 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 REAL 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 REAL 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.

Definition at line 85 of file ssterf.f.

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 *
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
real function slanst(NORM, N, D, E)
SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: slanst.f:100
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:63
subroutine slae2(A, B, C, RT1, RT2)
SLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition: slae2.f:102
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
Definition: slasrt.f:88
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function: