LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
slae2
subroutine slae2(a, b, c, rt1, rt2)
SLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition slae2.f:102
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:120
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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:100
slapy2
real function slapy2(x, y)
SLAPY2 returns sqrt(x2+y2).
Definition slapy2.f:63
slartg
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
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:143
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:110
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:199
slasrt
subroutine slasrt(id, n, d, info)
SLASRT sorts numbers in increasing or decreasing order.
Definition slasrt.f:88
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: