178 SUBROUTINE sstedc( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
187 INTEGER INFO, LDZ, LIWORK, LWORK, N
191 REAL D( * ), E( * ), WORK( * ), Z( LDZ, * )
198 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
202 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
203 $ lwmin, m, smlsiz, start, storez, strtrw
204 REAL EPS, ORGNRM, P, TINY
209 REAL SLAMCH, SLANST, SROUNDUP_LWORK
210 EXTERNAL ilaenv, lsame, slamch, slanst,
219 INTRINSIC abs, int, log, max, mod, real, sqrt
226 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
228 IF( lsame( compz,
'N' ) )
THEN
230 ELSE IF( lsame( compz,
'V' ) )
THEN
232 ELSE IF( lsame( compz,
'I' ) )
THEN
237 IF( icompz.LT.0 )
THEN
239 ELSE IF( n.LT.0 )
THEN
241 ELSE IF( ( ldz.LT.1 ) .OR.
242 $ ( icompz.GT.0 .AND. ldz.LT.max( 1, n ) ) )
THEN
250 smlsiz = ilaenv( 9,
'SSTEDC',
' ', 0, 0, 0, 0 )
251 IF( n.LE.1 .OR. icompz.EQ.0 )
THEN
254 ELSE IF( n.LE.smlsiz )
THEN
258 lgn = int( log( real( n ) )/log( two ) )
263 IF( icompz.EQ.1 )
THEN
264 lwmin = 1 + 3*n + 2*n*lgn + 4*n**2
265 liwmin = 6 + 6*n + 5*n*lgn
266 ELSE IF( icompz.EQ.2 )
THEN
267 lwmin = 1 + 4*n + n**2
271 work( 1 ) = sroundup_lwork(lwmin)
274 IF( lwork.LT.lwmin .AND. .NOT. lquery )
THEN
276 ELSE IF( liwork.LT.liwmin .AND. .NOT. lquery )
THEN
282 CALL xerbla(
'SSTEDC', -info )
284 ELSE IF (lquery)
THEN
309 IF( icompz.EQ.0 )
THEN
310 CALL ssterf( n, d, e, info )
317 IF( n.LE.smlsiz )
THEN
319 CALL ssteqr( compz, n, d, e, z, ldz, work, info )
326 IF( icompz.EQ.1 )
THEN
332 IF( icompz.EQ.2 )
THEN
333 CALL slaset(
'Full', n, n, zero, one, z, ldz )
338 orgnrm = slanst(
'M', n, d, e )
342 eps = slamch(
'Epsilon' )
349 IF( start.LE.n )
THEN
359 IF( finish.LT.n )
THEN
360 tiny = eps*sqrt( abs( d( finish ) ) )*
361 $ sqrt( abs( d( finish+1 ) ) )
362 IF( abs( e( finish ) ).GT.tiny )
THEN
370 m = finish - start + 1
375 IF( m.GT.smlsiz )
THEN
379 orgnrm = slanst(
'M', m, d( start ), e( start ) )
380 CALL slascl(
'G', 0, 0, orgnrm, one, m, 1, d( start ),
383 CALL slascl(
'G', 0, 0, orgnrm, one, m-1, 1,
387 IF( icompz.EQ.1 )
THEN
392 CALL slaed0( icompz, n, m, d( start ), e( start ),
393 $ z( strtrw, start ), ldz, work( 1 ), n,
394 $ work( storez ), iwork, info )
396 info = ( info / ( m+1 )+start-1 )*( n+1 ) +
397 $ mod( info, ( m+1 ) ) + start - 1
403 CALL slascl(
'G', 0, 0, one, orgnrm, m, 1, d( start ),
408 IF( icompz.EQ.1 )
THEN
414 CALL ssteqr(
'I', m, d( start ), e( start ), work,
416 $ work( m*m+1 ), info )
417 CALL slacpy(
'A', n, m, z( 1, start ), ldz,
418 $ work( storez ), n )
419 CALL sgemm(
'N',
'N', n, m, m, one,
420 $ work( storez ), n, work, m, zero,
421 $ z( 1, start ), ldz )
422 ELSE IF( icompz.EQ.2 )
THEN
423 CALL ssteqr(
'I', m, d( start ), e( start ),
424 $ z( start, start ), ldz, work, info )
426 CALL ssterf( m, d( start ), e( start ), info )
429 info = start*( n+1 ) + finish
440 IF( icompz.EQ.0 )
THEN
444 CALL slasrt(
'I', n, d, info )
455 IF( d( j ).LT.p )
THEN
463 CALL sswap( n, z( 1, i ), 1, z( 1, k ), 1 )
470 work( 1 ) = sroundup_lwork(lwmin)