186 SUBROUTINE sstedc( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
195 INTEGER INFO, LDZ, LIWORK, LWORK, N
199 REAL D( * ), E( * ), WORK( * ), Z( LDZ, * )
206 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
210 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
211 $ lwmin, m, smlsiz, start, storez, strtrw
212 REAL EPS, ORGNRM, P, TINY
218 EXTERNAL ilaenv, lsame, slamch, slanst
225 INTRINSIC abs, int, log, max, mod, real, sqrt
232 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
234 IF( lsame( compz,
'N' ) )
THEN
236 ELSE IF( lsame( compz,
'V' ) )
THEN
238 ELSE IF( lsame( compz,
'I' ) )
THEN
243 IF( icompz.LT.0 )
THEN
245 ELSE IF( n.LT.0 )
THEN
247 ELSE IF( ( ldz.LT.1 ) .OR.
248 $ ( icompz.GT.0 .AND. ldz.LT.max( 1, n ) ) )
THEN
256 smlsiz = ilaenv( 9,
'SSTEDC',
' ', 0, 0, 0, 0 )
257 IF( n.LE.1 .OR. icompz.EQ.0 )
THEN
260 ELSE IF( n.LE.smlsiz )
THEN
264 lgn = int( log( real( n ) )/log( two ) )
269 IF( icompz.EQ.1 )
THEN
270 lwmin = 1 + 3*n + 2*n*lgn + 4*n**2
271 liwmin = 6 + 6*n + 5*n*lgn
272 ELSE IF( icompz.EQ.2 )
THEN
273 lwmin = 1 + 4*n + n**2
280 IF( lwork.LT.lwmin .AND. .NOT. lquery )
THEN
282 ELSE IF( liwork.LT.liwmin .AND. .NOT. lquery )
THEN
288 CALL xerbla(
'SSTEDC', -info )
290 ELSE IF (lquery)
THEN
315 IF( icompz.EQ.0 )
THEN
316 CALL ssterf( n, d, e, info )
323 IF( n.LE.smlsiz )
THEN
325 CALL ssteqr( compz, n, d, e, z, ldz, work, info )
332 IF( icompz.EQ.1 )
THEN
338 IF( icompz.EQ.2 )
THEN
339 CALL slaset(
'Full', n, n, zero, one, z, ldz )
344 orgnrm = slanst(
'M', n, d, e )
348 eps = slamch(
'Epsilon' )
355 IF( start.LE.n )
THEN
365 IF( finish.LT.n )
THEN
366 tiny = eps*sqrt( abs( d( finish ) ) )*
367 $ sqrt( abs( d( finish+1 ) ) )
368 IF( abs( e( finish ) ).GT.tiny )
THEN
376 m = finish - start + 1
381 IF( m.GT.smlsiz )
THEN
385 orgnrm = slanst(
'M', m, d( start ), e( start ) )
386 CALL slascl(
'G', 0, 0, orgnrm, one, m, 1, d( start ), m,
388 CALL slascl(
'G', 0, 0, orgnrm, one, m-1, 1, e( start ),
391 IF( icompz.EQ.1 )
THEN
396 CALL slaed0( icompz, n, m, d( start ), e( start ),
397 $ z( strtrw, start ), ldz, work( 1 ), n,
398 $ work( storez ), iwork, info )
400 info = ( info / ( m+1 )+start-1 )*( n+1 ) +
401 $ mod( info, ( m+1 ) ) + start - 1
407 CALL slascl(
'G', 0, 0, one, orgnrm, m, 1, d( start ), m,
411 IF( icompz.EQ.1 )
THEN
417 CALL ssteqr(
'I', m, d( start ), e( start ), work, m,
418 $ work( m*m+1 ), info )
419 CALL slacpy(
'A', n, m, z( 1, start ), ldz,
420 $ work( storez ), n )
421 CALL sgemm(
'N',
'N', n, m, m, one,
422 $ work( storez ), n, work, m, zero,
423 $ z( 1, start ), ldz )
424 ELSE IF( icompz.EQ.2 )
THEN
425 CALL ssteqr(
'I', m, d( start ), e( start ),
426 $ z( start, start ), ldz, work, info )
428 CALL ssterf( m, d( start ), e( start ), info )
431 info = start*( n+1 ) + finish
442 IF( icompz.EQ.0 )
THEN
446 CALL slasrt(
'I', n, d, info )
457 IF( d( j ).LT.p )
THEN
465 CALL sswap( n, z( 1, i ), 1, z( 1, k ), 1 )
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.
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.
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
subroutine sstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEDC
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
subroutine ssterf(N, D, E, INFO)
SSTERF
subroutine slaed0(ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, WORK, IWORK, INFO)
SLAED0 used by SSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM