LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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 steqr
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,
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 xerbla(srname, info)
Definition cblat2.f:3285
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 slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
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 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 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 slasrt(id, n, d, info)
SLASRT sorts numbers in increasing or decreasing order.
Definition slasrt.f:88
subroutine ssteqr(compz, n, d, e, z, ldz, work, info)
SSTEQR
Definition ssteqr.f:131
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82