LAPACK  3.10.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.

Definition at line 130 of file ssteqr.f.

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