174 SUBROUTINE sstein( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
175 $ iwork, ifail, info )
183 INTEGER INFO, LDZ, M, N
186 INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
188 REAL D( * ), E( * ), W( * ), WORK( * ), Z( ldz, * )
194 REAL ZERO, ONE, TEN, ODM3, ODM1
195 parameter ( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+1,
196 $ odm3 = 1.0e-3, odm1 = 1.0e-1 )
197 INTEGER MAXITS, EXTRA
198 parameter ( maxits = 5, extra = 2 )
201 INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
202 $ indrv2, indrv3, indrv4, indrv5, its, j, j1,
203 $ jblk, jmax, nblk, nrmchk
204 REAL CTR, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
205 $ scl, sep, stpcrt, tol, xj, xjm
212 REAL SASUM, SDOT, SLAMCH, SNRM2
213 EXTERNAL isamax, sasum, sdot, slamch, snrm2
220 INTRINSIC abs, max, sqrt
233 ELSE IF( m.LT.0 .OR. m.GT.n )
THEN
235 ELSE IF( ldz.LT.max( 1, n ) )
THEN
239 IF( iblock( j ).LT.iblock( j-1 ) )
THEN
243 IF( iblock( j ).EQ.iblock( j-1 ) .AND. w( j ).LT.w( j-1 ) )
253 CALL xerbla(
'SSTEIN', -info )
259 IF( n.EQ.0 .OR. m.EQ.0 )
THEN
261 ELSE IF( n.EQ.1 )
THEN
268 eps = slamch(
'Precision' )
287 DO 160 nblk = 1, iblock( m )
294 b1 = isplit( nblk-1 ) + 1
304 onenrm = abs( d( b1 ) ) + abs( e( b1 ) )
305 onenrm = max( onenrm, abs( d( bn ) )+abs( e( bn-1 ) ) )
306 DO 50 i = b1 + 1, bn - 1
307 onenrm = max( onenrm, abs( d( i ) )+abs( e( i-1 ) )+
312 stpcrt = sqrt( odm1 / blksiz )
319 IF( iblock( j ).NE.nblk )
THEN
328 IF( blksiz.EQ.1 )
THEN
329 work( indrv1+1 ) = one
349 CALL slarnv( 2, iseed, blksiz, work( indrv1+1 ) )
353 CALL scopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
354 CALL scopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
355 CALL scopy( blksiz-1, e( b1 ), 1, work( indrv3+1 ), 1 )
360 CALL slagtf( blksiz, work( indrv4+1 ), xj, work( indrv2+2 ),
361 $ work( indrv3+1 ), tol, work( indrv5+1 ), iwork,
373 jmax = isamax( blksiz, work( indrv1+1 ), 1 )
374 scl = blksiz*onenrm*max( eps,
375 $ abs( work( indrv4+blksiz ) ) ) /
376 $ abs( work( indrv1+jmax ) )
377 CALL sscal( blksiz, scl, work( indrv1+1 ), 1 )
381 CALL slagts( -1, blksiz, work( indrv4+1 ), work( indrv2+2 ),
382 $ work( indrv3+1 ), work( indrv5+1 ), iwork,
383 $ work( indrv1+1 ), tol, iinfo )
390 IF( abs( xj-xjm ).GT.ortol )
392 IF( gpind.NE.j )
THEN
393 DO 80 i = gpind, j - 1
394 ctr = -sdot( blksiz, work( indrv1+1 ), 1, z( b1, i ),
396 CALL saxpy( blksiz, ctr, z( b1, i ), 1,
397 $ work( indrv1+1 ), 1 )
404 jmax = isamax( blksiz, work( indrv1+1 ), 1 )
405 nrm = abs( work( indrv1+jmax ) )
413 IF( nrmchk.LT.extra+1 )
428 scl = one / snrm2( blksiz, work( indrv1+1 ), 1 )
429 jmax = isamax( blksiz, work( indrv1+1 ), 1 )
430 IF( work( indrv1+jmax ).LT.zero )
432 CALL sscal( blksiz, scl, work( indrv1+1 ), 1 )
438 z( b1+i-1, j ) = work( indrv1+i )
subroutine sstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
SSTEIN
subroutine slarnv(IDIST, ISEED, N, X)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine slagtf(N, A, LAMBDA, B, C, TOL, D, IN, INFO)
SLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix...
subroutine slagts(JOB, N, A, B, C, D, IN, Y, TOL, INFO)
SLAGTS solves the system of equations (T-λI)x = y or (T-λI)Tx = y,where T is a general tridiagonal ma...
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY