LAPACK  3.9.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.
Date
December 2016

Definition at line 88 of file ssterf.f.

88 *
89 * -- LAPACK computational routine (version 3.7.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 * December 2016
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 *
Here is the call graph for this function:
slae2
subroutine slae2(A, B, C, RT1, RT2)
SLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition: slae2.f:104
slanst
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:102
slascl
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:145
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
slapy2
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:65
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:70
slasrt
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
Definition: slasrt.f:90