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

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

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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