LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
ssteqr.f
Go to the documentation of this file.
1 *> \brief \b SSTEQR
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SSTEQR + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssteqr.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssteqr.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssteqr.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER COMPZ
25 * INTEGER INFO, LDZ, N
26 * ..
27 * .. Array Arguments ..
28 * REAL D( * ), E( * ), WORK( * ), Z( LDZ, * )
29 * ..
30 *
31 *
32 *> \par Purpose:
33 * =============
34 *>
35 *> \verbatim
36 *>
37 *> SSTEQR computes all eigenvalues and, optionally, eigenvectors of a
38 *> symmetric tridiagonal matrix using the implicit QL or QR method.
39 *> The eigenvectors of a full or band symmetric matrix can also be found
40 *> if SSYTRD or SSPTRD or SSBTRD has been used to reduce this matrix to
41 *> tridiagonal form.
42 *> \endverbatim
43 *
44 * Arguments:
45 * ==========
46 *
47 *> \param[in] COMPZ
48 *> \verbatim
49 *> COMPZ is CHARACTER*1
50 *> = 'N': Compute eigenvalues only.
51 *> = 'V': Compute eigenvalues and eigenvectors of the original
52 *> symmetric matrix. On entry, Z must contain the
53 *> orthogonal matrix used to reduce the original matrix
54 *> to tridiagonal form.
55 *> = 'I': Compute eigenvalues and eigenvectors of the
56 *> tridiagonal matrix. Z is initialized to the identity
57 *> matrix.
58 *> \endverbatim
59 *>
60 *> \param[in] N
61 *> \verbatim
62 *> N is INTEGER
63 *> The order of the matrix. N >= 0.
64 *> \endverbatim
65 *>
66 *> \param[in,out] D
67 *> \verbatim
68 *> D is REAL array, dimension (N)
69 *> On entry, the diagonal elements of the tridiagonal matrix.
70 *> On exit, if INFO = 0, the eigenvalues in ascending order.
71 *> \endverbatim
72 *>
73 *> \param[in,out] E
74 *> \verbatim
75 *> E is REAL array, dimension (N-1)
76 *> On entry, the (n-1) subdiagonal elements of the tridiagonal
77 *> matrix.
78 *> On exit, E has been destroyed.
79 *> \endverbatim
80 *>
81 *> \param[in,out] Z
82 *> \verbatim
83 *> Z is REAL array, dimension (LDZ, N)
84 *> On entry, if COMPZ = 'V', then Z contains the orthogonal
85 *> matrix used in the reduction to tridiagonal form.
86 *> On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
87 *> orthonormal eigenvectors of the original symmetric matrix,
88 *> and if COMPZ = 'I', Z contains the orthonormal eigenvectors
89 *> of the symmetric tridiagonal matrix.
90 *> If COMPZ = 'N', then Z is not referenced.
91 *> \endverbatim
92 *>
93 *> \param[in] LDZ
94 *> \verbatim
95 *> LDZ is INTEGER
96 *> The leading dimension of the array Z. LDZ >= 1, and if
97 *> eigenvectors are desired, then LDZ >= max(1,N).
98 *> \endverbatim
99 *>
100 *> \param[out] WORK
101 *> \verbatim
102 *> WORK is REAL array, dimension (max(1,2*N-2))
103 *> If COMPZ = 'N', then WORK is not referenced.
104 *> \endverbatim
105 *>
106 *> \param[out] INFO
107 *> \verbatim
108 *> INFO is INTEGER
109 *> = 0: successful exit
110 *> < 0: if INFO = -i, the i-th argument had an illegal value
111 *> > 0: the algorithm has failed to find all the eigenvalues in
112 *> a total of 30*N iterations; if INFO = i, then i
113 *> elements of E have not converged to zero; on exit, D
114 *> and E contain the elements of a symmetric tridiagonal
115 *> matrix which is orthogonally similar to the original
116 *> matrix.
117 *> \endverbatim
118 *
119 * Authors:
120 * ========
121 *
122 *> \author Univ. of Tennessee
123 *> \author Univ. of California Berkeley
124 *> \author Univ. of Colorado Denver
125 *> \author NAG Ltd.
126 *
127 *> \ingroup auxOTHERcomputational
128 *
129 * =====================================================================
130  SUBROUTINE ssteqr( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
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 *
569  END
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
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
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
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:131
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