188 SUBROUTINE sstedc( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
198 INTEGER INFO, LDZ, LIWORK, LWORK, N
202 REAL D( * ), E( * ), WORK( * ), Z( ldz, * )
209 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
213 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
214 $ lwmin, m, smlsiz, start, storez, strtrw
215 REAL EPS, ORGNRM, P, TINY
221 EXTERNAL ilaenv, lsame, slamch, slanst
228 INTRINSIC abs, int, log, max, mod,
REAL, SQRT
235 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
237 IF( lsame( compz,
'N' ) )
THEN
239 ELSE IF( lsame( compz,
'V' ) )
THEN
241 ELSE IF( lsame( compz,
'I' ) )
THEN
246 IF( icompz.LT.0 )
THEN
248 ELSE IF( n.LT.0 )
THEN
250 ELSE IF( ( ldz.LT.1 ) .OR.
251 $ ( icompz.GT.0 .AND. ldz.LT.max( 1, n ) ) )
THEN
259 smlsiz = ilaenv( 9,
'SSTEDC',
' ', 0, 0, 0, 0 )
260 IF( n.LE.1 .OR. icompz.EQ.0 )
THEN
263 ELSE IF( n.LE.smlsiz )
THEN
267 lgn = int( log(
REAL( N ) )/log( TWO ) )
272 IF( icompz.EQ.1 )
THEN
273 lwmin = 1 + 3*n + 2*n*lgn + 4*n**2
274 liwmin = 6 + 6*n + 5*n*lgn
275 ELSE IF( icompz.EQ.2 )
THEN
276 lwmin = 1 + 4*n + n**2
283 IF( lwork.LT.lwmin .AND. .NOT. lquery )
THEN
285 ELSE IF( liwork.LT.liwmin .AND. .NOT. lquery )
THEN
291 CALL
xerbla(
'SSTEDC', -info )
293 ELSE IF (lquery)
THEN
318 IF( icompz.EQ.0 )
THEN
319 CALL
ssterf( n, d, e, info )
326 IF( n.LE.smlsiz )
THEN
328 CALL
ssteqr( compz, n, d, e, z, ldz, work, info )
335 IF( icompz.EQ.1 )
THEN
341 IF( icompz.EQ.2 )
THEN
342 CALL
slaset(
'Full', n, n, zero, one, z, ldz )
347 orgnrm = slanst(
'M', n, d, e )
351 eps = slamch(
'Epsilon' )
358 IF( start.LE.n )
THEN
368 IF( finish.LT.n )
THEN
369 tiny = eps*sqrt( abs( d( finish ) ) )*
370 $ sqrt( abs( d( finish+1 ) ) )
371 IF( abs( e( finish ) ).GT.tiny )
THEN
379 m = finish - start + 1
384 IF( m.GT.smlsiz )
THEN
388 orgnrm = slanst(
'M', m, d( start ), e( start ) )
389 CALL
slascl(
'G', 0, 0, orgnrm, one, m, 1, d( start ), m,
391 CALL
slascl(
'G', 0, 0, orgnrm, one, m-1, 1, e( start ),
394 IF( icompz.EQ.1 )
THEN
399 CALL
slaed0( icompz, n, m, d( start ), e( start ),
400 $ z( strtrw, start ), ldz, work( 1 ), n,
401 $ work( storez ), iwork, info )
403 info = ( info / ( m+1 )+start-1 )*( n+1 ) +
404 $ mod( info, ( m+1 ) ) + start - 1
410 CALL
slascl(
'G', 0, 0, one, orgnrm, m, 1, d( start ), m,
414 IF( icompz.EQ.1 )
THEN
420 CALL
ssteqr(
'I', m, d( start ), e( start ), work, m,
421 $ work( m*m+1 ), info )
422 CALL
slacpy(
'A', n, m, z( 1, start ), ldz,
423 $ work( storez ), n )
424 CALL
sgemm(
'N',
'N', n, m, m, one,
425 $ work( storez ), n, work, m, zero,
426 $ z( 1, start ), ldz )
427 ELSE IF( icompz.EQ.2 )
THEN
428 CALL
ssteqr(
'I', m, d( start ), e( start ),
429 $ z( start, start ), ldz, work, info )
431 CALL
ssterf( m, d( start ), e( start ), info )
434 info = start*( n+1 ) + finish
450 IF( icompz.EQ.0 )
THEN
454 CALL
slasrt(
'I', n, d, info )
465 IF( d( j ).LT.p )
THEN
473 CALL
sswap( n, z( 1, i ), 1, z( 1, k ), 1 )
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 sswap(N, SX, INCX, SY, INCY)
SSWAP
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 xerbla(SRNAME, INFO)
XERBLA
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEBZ
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
subroutine ssterf(N, D, E, INFO)
SSTERF
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR