 LAPACK  3.9.0 LAPACK: Linear Algebra PACKage

## ◆ ssteqr()

 subroutine ssteqr ( character COMPZ, integer N, real, dimension( * ) D, real, dimension( * ) E, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer INFO )

SSTEQR

Purpose:
``` SSTEQR computes all eigenvalues and, optionally, eigenvectors of a
symmetric tridiagonal matrix using the implicit QL or QR method.
The eigenvectors of a full or band symmetric matrix can also be found
if SSYTRD or SSPTRD or SSBTRD has been used to reduce this matrix to
tridiagonal form.```
Parameters
 [in] COMPZ ``` COMPZ is CHARACTER*1 = 'N': Compute eigenvalues only. = 'V': Compute eigenvalues and eigenvectors of the original symmetric matrix. On entry, Z must contain the orthogonal matrix used to reduce the original matrix to tridiagonal form. = 'I': Compute eigenvalues and eigenvectors of the tridiagonal matrix. Z is initialized to the identity matrix.``` [in] N ``` N is INTEGER The order of the matrix. N >= 0.``` [in,out] D ``` D is REAL array, dimension (N) On entry, the 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.``` [in,out] Z ``` Z is REAL array, dimension (LDZ, N) On entry, if COMPZ = 'V', then Z contains the orthogonal matrix used in the reduction to tridiagonal form. On exit, if INFO = 0, then if COMPZ = 'V', Z contains the orthonormal eigenvectors of the original symmetric matrix, and if COMPZ = 'I', Z contains the orthonormal eigenvectors of the symmetric tridiagonal matrix. If COMPZ = 'N', then Z is not referenced.``` [in] LDZ ``` LDZ is INTEGER The leading dimension of the array Z. LDZ >= 1, and if eigenvectors are desired, then LDZ >= max(1,N).``` [out] WORK ``` WORK is REAL array, dimension (max(1,2*N-2)) If COMPZ = 'N', then WORK is not referenced.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: the algorithm has failed to find all the eigenvalues in a total of 30*N iterations; if INFO = i, then i elements of E have not converged to zero; on exit, D and E contain the elements of a symmetric tridiagonal matrix which is orthogonally similar to the original matrix.```
Date
December 2016

Definition at line 133 of file ssteqr.f.

133 *
134 * -- LAPACK computational routine (version 3.7.0) --
135 * -- LAPACK is a software package provided by Univ. of Tennessee, --
136 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
137 * December 2016
138 *
139 * .. Scalar Arguments ..
140  CHARACTER COMPZ
141  INTEGER INFO, LDZ, N
142 * ..
143 * .. Array Arguments ..
144  REAL D( * ), E( * ), WORK( * ), Z( LDZ, * )
145 * ..
146 *
147 * =====================================================================
148 *
149 * .. Parameters ..
150  REAL ZERO, ONE, TWO, THREE
151  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
152  \$ three = 3.0e0 )
153  INTEGER MAXIT
154  parameter( maxit = 30 )
155 * ..
156 * .. Local Scalars ..
157  INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
158  \$ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
159  \$ NM1, NMAXIT
160  REAL ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
161  \$ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
162 * ..
163 * .. External Functions ..
164  LOGICAL LSAME
165  REAL SLAMCH, SLANST, SLAPY2
166  EXTERNAL lsame, slamch, slanst, slapy2
167 * ..
168 * .. External Subroutines ..
169  EXTERNAL slae2, slaev2, slartg, slascl, slaset, slasr,
170  \$ slasrt, sswap, xerbla
171 * ..
172 * .. Intrinsic Functions ..
173  INTRINSIC abs, max, sign, sqrt
174 * ..
175 * .. Executable Statements ..
176 *
177 * Test the input parameters.
178 *
179  info = 0
180 *
181  IF( lsame( compz, 'N' ) ) THEN
182  icompz = 0
183  ELSE IF( lsame( compz, 'V' ) ) THEN
184  icompz = 1
185  ELSE IF( lsame( compz, 'I' ) ) THEN
186  icompz = 2
187  ELSE
188  icompz = -1
189  END IF
190  IF( icompz.LT.0 ) THEN
191  info = -1
192  ELSE IF( n.LT.0 ) THEN
193  info = -2
194  ELSE IF( ( ldz.LT.1 ) .OR. ( icompz.GT.0 .AND. ldz.LT.max( 1,
195  \$ n ) ) ) THEN
196  info = -6
197  END IF
198  IF( info.NE.0 ) THEN
199  CALL xerbla( 'SSTEQR', -info )
200  RETURN
201  END IF
202 *
203 * Quick return if possible
204 *
205  IF( n.EQ.0 )
206  \$ RETURN
207 *
208  IF( n.EQ.1 ) THEN
209  IF( icompz.EQ.2 )
210  \$ z( 1, 1 ) = one
211  RETURN
212  END IF
213 *
214 * Determine the unit roundoff and over/underflow thresholds.
215 *
216  eps = slamch( 'E' )
217  eps2 = eps**2
218  safmin = slamch( 'S' )
219  safmax = one / safmin
220  ssfmax = sqrt( safmax ) / three
221  ssfmin = sqrt( safmin ) / eps2
222 *
223 * Compute the eigenvalues and eigenvectors of the tridiagonal
224 * matrix.
225 *
226  IF( icompz.EQ.2 )
227  \$ CALL slaset( 'Full', n, n, zero, one, z, ldz )
228 *
229  nmaxit = n*maxit
230  jtot = 0
231 *
232 * Determine where the matrix splits and choose QL or QR iteration
233 * for each block, according to whether top or bottom diagonal
234 * element is smaller.
235 *
236  l1 = 1
237  nm1 = n - 1
238 *
239  10 CONTINUE
240  IF( l1.GT.n )
241  \$ GO TO 160
242  IF( l1.GT.1 )
243  \$ e( l1-1 ) = zero
244  IF( l1.LE.nm1 ) THEN
245  DO 20 m = l1, nm1
246  tst = abs( e( m ) )
247  IF( tst.EQ.zero )
248  \$ GO TO 30
249  IF( tst.LE.( sqrt( abs( d( m ) ) )*sqrt( abs( d( m+
250  \$ 1 ) ) ) )*eps ) THEN
251  e( m ) = zero
252  GO TO 30
253  END IF
254  20 CONTINUE
255  END IF
256  m = n
257 *
258  30 CONTINUE
259  l = l1
260  lsv = l
261  lend = m
262  lendsv = lend
263  l1 = m + 1
264  IF( lend.EQ.l )
265  \$ GO TO 10
266 *
267 * Scale submatrix in rows and columns L to LEND
268 *
269  anorm = slanst( 'M', lend-l+1, d( l ), e( l ) )
270  iscale = 0
271  IF( anorm.EQ.zero )
272  \$ GO TO 10
273  IF( anorm.GT.ssfmax ) THEN
274  iscale = 1
275  CALL slascl( 'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
276  \$ info )
277  CALL slascl( 'G', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
278  \$ info )
279  ELSE IF( anorm.LT.ssfmin ) THEN
280  iscale = 2
281  CALL slascl( 'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
282  \$ info )
283  CALL slascl( 'G', 0, 0, anorm, ssfmin, lend-l, 1, e( l ), n,
284  \$ info )
285  END IF
286 *
287 * Choose between QL and QR iteration
288 *
289  IF( abs( d( lend ) ).LT.abs( d( l ) ) ) THEN
290  lend = lsv
291  l = lendsv
292  END IF
293 *
294  IF( lend.GT.l ) THEN
295 *
296 * QL Iteration
297 *
298 * Look for small subdiagonal element.
299 *
300  40 CONTINUE
301  IF( l.NE.lend ) THEN
302  lendm1 = lend - 1
303  DO 50 m = l, lendm1
304  tst = abs( e( m ) )**2
305  IF( tst.LE.( eps2*abs( d( m ) ) )*abs( d( m+1 ) )+
306  \$ safmin )GO TO 60
307  50 CONTINUE
308  END IF
309 *
310  m = lend
311 *
312  60 CONTINUE
313  IF( m.LT.lend )
314  \$ e( m ) = zero
315  p = d( l )
316  IF( m.EQ.l )
317  \$ GO TO 80
318 *
319 * If remaining matrix is 2-by-2, use SLAE2 or SLAEV2
320 * to compute its eigensystem.
321 *
322  IF( m.EQ.l+1 ) THEN
323  IF( icompz.GT.0 ) THEN
324  CALL slaev2( d( l ), e( l ), d( l+1 ), rt1, rt2, c, s )
325  work( l ) = c
326  work( n-1+l ) = s
327  CALL slasr( 'R', 'V', 'B', n, 2, work( l ),
328  \$ work( n-1+l ), z( 1, l ), ldz )
329  ELSE
330  CALL slae2( d( l ), e( l ), d( l+1 ), rt1, rt2 )
331  END IF
332  d( l ) = rt1
333  d( l+1 ) = rt2
334  e( l ) = zero
335  l = l + 2
336  IF( l.LE.lend )
337  \$ GO TO 40
338  GO TO 140
339  END IF
340 *
341  IF( jtot.EQ.nmaxit )
342  \$ GO TO 140
343  jtot = jtot + 1
344 *
345 * Form shift.
346 *
347  g = ( d( l+1 )-p ) / ( two*e( l ) )
348  r = slapy2( g, one )
349  g = d( m ) - p + ( e( l ) / ( g+sign( r, g ) ) )
350 *
351  s = one
352  c = one
353  p = zero
354 *
355 * Inner loop
356 *
357  mm1 = m - 1
358  DO 70 i = mm1, l, -1
359  f = s*e( i )
360  b = c*e( i )
361  CALL slartg( g, f, c, s, r )
362  IF( i.NE.m-1 )
363  \$ e( i+1 ) = r
364  g = d( i+1 ) - p
365  r = ( d( i )-g )*s + two*c*b
366  p = s*r
367  d( i+1 ) = g + p
368  g = c*r - b
369 *
370 * If eigenvectors are desired, then save rotations.
371 *
372  IF( icompz.GT.0 ) THEN
373  work( i ) = c
374  work( n-1+i ) = -s
375  END IF
376 *
377  70 CONTINUE
378 *
379 * If eigenvectors are desired, then apply saved rotations.
380 *
381  IF( icompz.GT.0 ) THEN
382  mm = m - l + 1
383  CALL slasr( 'R', 'V', 'B', n, mm, work( l ), work( n-1+l ),
384  \$ z( 1, l ), ldz )
385  END IF
386 *
387  d( l ) = d( l ) - p
388  e( l ) = g
389  GO TO 40
390 *
391 * Eigenvalue found.
392 *
393  80 CONTINUE
394  d( l ) = p
395 *
396  l = l + 1
397  IF( l.LE.lend )
398  \$ GO TO 40
399  GO TO 140
400 *
401  ELSE
402 *
403 * QR Iteration
404 *
405 * Look for small superdiagonal element.
406 *
407  90 CONTINUE
408  IF( l.NE.lend ) THEN
409  lendp1 = lend + 1
410  DO 100 m = l, lendp1, -1
411  tst = abs( e( m-1 ) )**2
412  IF( tst.LE.( eps2*abs( d( m ) ) )*abs( d( m-1 ) )+
413  \$ safmin )GO TO 110
414  100 CONTINUE
415  END IF
416 *
417  m = lend
418 *
419  110 CONTINUE
420  IF( m.GT.lend )
421  \$ e( m-1 ) = zero
422  p = d( l )
423  IF( m.EQ.l )
424  \$ GO TO 130
425 *
426 * If remaining matrix is 2-by-2, use SLAE2 or SLAEV2
427 * to compute its eigensystem.
428 *
429  IF( m.EQ.l-1 ) THEN
430  IF( icompz.GT.0 ) THEN
431  CALL slaev2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2, c, s )
432  work( m ) = c
433  work( n-1+m ) = s
434  CALL slasr( 'R', 'V', 'F', n, 2, work( m ),
435  \$ work( n-1+m ), z( 1, l-1 ), ldz )
436  ELSE
437  CALL slae2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2 )
438  END IF
439  d( l-1 ) = rt1
440  d( l ) = rt2
441  e( l-1 ) = zero
442  l = l - 2
443  IF( l.GE.lend )
444  \$ GO TO 90
445  GO TO 140
446  END IF
447 *
448  IF( jtot.EQ.nmaxit )
449  \$ GO TO 140
450  jtot = jtot + 1
451 *
452 * Form shift.
453 *
454  g = ( d( l-1 )-p ) / ( two*e( l-1 ) )
455  r = slapy2( g, one )
456  g = d( m ) - p + ( e( l-1 ) / ( g+sign( r, g ) ) )
457 *
458  s = one
459  c = one
460  p = zero
461 *
462 * Inner loop
463 *
464  lm1 = l - 1
465  DO 120 i = m, lm1
466  f = s*e( i )
467  b = c*e( i )
468  CALL slartg( g, f, c, s, r )
469  IF( i.NE.m )
470  \$ e( i-1 ) = r
471  g = d( i ) - p
472  r = ( d( i+1 )-g )*s + two*c*b
473  p = s*r
474  d( i ) = g + p
475  g = c*r - b
476 *
477 * If eigenvectors are desired, then save rotations.
478 *
479  IF( icompz.GT.0 ) THEN
480  work( i ) = c
481  work( n-1+i ) = s
482  END IF
483 *
484  120 CONTINUE
485 *
486 * If eigenvectors are desired, then apply saved rotations.
487 *
488  IF( icompz.GT.0 ) THEN
489  mm = l - m + 1
490  CALL slasr( 'R', 'V', 'F', n, mm, work( m ), work( n-1+m ),
491  \$ z( 1, m ), ldz )
492  END IF
493 *
494  d( l ) = d( l ) - p
495  e( lm1 ) = g
496  GO TO 90
497 *
498 * Eigenvalue found.
499 *
500  130 CONTINUE
501  d( l ) = p
502 *
503  l = l - 1
504  IF( l.GE.lend )
505  \$ GO TO 90
506  GO TO 140
507 *
508  END IF
509 *
510 * Undo scaling if necessary
511 *
512  140 CONTINUE
513  IF( iscale.EQ.1 ) THEN
514  CALL slascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
515  \$ d( lsv ), n, info )
516  CALL slascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e( lsv ),
517  \$ n, info )
518  ELSE IF( iscale.EQ.2 ) THEN
519  CALL slascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
520  \$ d( lsv ), n, info )
521  CALL slascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv, 1, e( lsv ),
522  \$ n, info )
523  END IF
524 *
525 * Check for no convergence to an eigenvalue after a total
526 * of N*MAXIT iterations.
527 *
528  IF( jtot.LT.nmaxit )
529  \$ GO TO 10
530  DO 150 i = 1, n - 1
531  IF( e( i ).NE.zero )
532  \$ info = info + 1
533  150 CONTINUE
534  GO TO 190
535 *
536 * Order eigenvalues and eigenvectors.
537 *
538  160 CONTINUE
539  IF( icompz.EQ.0 ) THEN
540 *
541 * Use Quick Sort
542 *
543  CALL slasrt( 'I', n, d, info )
544 *
545  ELSE
546 *
547 * Use Selection Sort to minimize swaps of eigenvectors
548 *
549  DO 180 ii = 2, n
550  i = ii - 1
551  k = i
552  p = d( i )
553  DO 170 j = ii, n
554  IF( d( j ).LT.p ) THEN
555  k = j
556  p = d( j )
557  END IF
558  170 CONTINUE
559  IF( k.NE.i ) THEN
560  d( k ) = d( i )
561  d( i ) = p
562  CALL sswap( n, z( 1, i ), 1, z( 1, k ), 1 )
563  END IF
564  180 CONTINUE
565  END IF
566 *
567  190 CONTINUE
568  RETURN
569 *
570 * End of SSTEQR
571 *
Here is the call graph for this function:
Here is the caller 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
slaev2
subroutine slaev2(A, B, C, RT1, RT2, CS1, SN1)
SLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.
Definition: slaev2.f:122
sswap
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:84
slasr
subroutine slasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
SLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition: slasr.f:201
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
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
slapy2
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:65
slaset
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:112
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:70
slartg
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99
slasrt
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
Definition: slasrt.f:90