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
445 IF( icompz.EQ.0 )
THEN
449 CALL slasrt(
'I', n, d, info )
460 IF( d( j ).LT.p )
THEN
468 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 sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
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 ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
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 ssterf(N, D, E, INFO)
SSTERF
subroutine sstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEDC